#include "kernel/driver.hpp"
#include "yocto/code/utils.hpp"

void Driver:: Init( Driver::Type &odeint, const Cell &cell )
{
	odeint.start( cell.nvar );
	odeint.eps  = cell.chemsys.ftol;
	odeint.hmin = 0;
	cell.dt_ctl = cell.dt_max * 1e-3;
}

void Driver:: Step( Driver::Type &odeint, Cell &cell, Real t_old, Real t_new )
{
	
	//--------------------------------------------------------------------------
	// Transfert S_in into variables, as well as Em
	//--------------------------------------------------------------------------
	Solution    &S_in  = cell.solutions[ Cell:: iS_in ];
	const size_t M     = cell.M;
	const size_t jEm   = M + Cell:: iEm;
	array<Real> &F     = cell.variables;
	
	S_in.save( F );
	F[ jEm ] = cell.Em;

	const Real dt     = cell.dt_max;
	Real       t_curr = t_old;
	cell.dt_ctl       = min_of(cell.dt_ctl,cell.dt_max);
	do 
	{
		//----------------------------------------------------------------------
		//-- compute next time
		//----------------------------------------------------------------------
		const Real t_next = max_of(t_curr+dt,t_new);
		
		//----------------------------------------------------------------------
		//-- step to t_next, adaptive
		//----------------------------------------------------------------------
		odeint( cell.drvs, F, t_curr, t_next, cell.dt_ctl );
		
		//----------------------------------------------------------------------
		//-- feed the chemical system
		//----------------------------------------------------------------------
		for( size_t i=M;i>0;--i )
		{
			Real &c = F[i];
			if( c <= 0 ) c = 0;
		}
		cell.chemsys.load_data( F );
		
		//----------------------------------------------------------------------
		//-- normalize it
		//----------------------------------------------------------------------
		cell.chemsys.solveGamma(t_next);
		
		//----------------------------------------------------------------------
		//-- store in current solution, update Em as well
		//----------------------------------------------------------------------
		S_in.load_from( cell.chemsys );
		cell.Em = F[ jEm ];
		t_curr = t_next;
	} while (t_curr<t_new);
	
    //--------------------------------------------------------------------------
    //
    //--------------------------------------------------------------------------
    
}

#if 0

//--------------------------------------------------------------------------
// Take a step
//--------------------------------------------------------------------------
odeint( cell.drvs, F, t_old, t_new, cell.dt_ctl );

//--------------------------------------------------------------------------
// Transfert F into chemsys
//--------------------------------------------------------------------------
cell.chemsys.load_data( F );

//--------------------------------------------------------------------------
// Newton Step to correct it
//--------------------------------------------------------------------------
cell.chemsys.solveGamma(t_new);

//--------------------------------------------------------------------------
// Transfer chemsys into S_in, as well as Em
//--------------------------------------------------------------------------
S_in.load_from( cell.chemsys );
cell.Em = F[ jEm ];
#endif
